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The Formation of Submillimetre-Bright Galaxies 
from Gas Infall over a Billion Years 

Desika Narayanan^Matthew Turk^, Robert Feldmann^’^^, Thomas Robitaille'^, Philip Hopkins^, Robert Thompson^’^, Christo¬ 
pher Hayward^’^, David Ball^’^, Claude-Andre Faucher-Giguere^^, Dusan Keres^^ 


Sulflnillimetre-luminous galaxies at high-redshift are the most lu¬ 
minous, heavily star-forming galaxies in the Univers^, and are 
characterised by prodigious emission in the far-infrared at 850 mi¬ 
crons {Ssbo > 5 mjy). They reside in halos ~ have 

low gas fractions compared to main sequence disks at a compara¬ 
ble redshifP, trace complex environmentJ3^ and are not easily ob¬ 
servable at optical wavelengthJ^ Their physical origin remains un¬ 
clear. Simulations have been able to form galaxies with the requi¬ 
site luminosities, but have otherwise been unable to simultaneously 
match the st ellar masses, star formation rates, gas fractions and 
environment^®^. Here we report a cosmological hydrodynamic 
galaxy formation simulation that is able to form a submillimetre 
galaxy which simultaneously satisfies the broad range of observed 
physical constraints. We find that groups of galaxies residing in 
massive dark matter halos have rising star formation histories that 
peak at collective rates ~ 500 — 10^Moyr“^ at z = 2 — 3, by which 
time the interstellar medium is sufficiently enriched with metals 
that the region may be observed as a submillimetre-selected sys¬ 
tem. The intense star formation rates are fueled in part by a reser¬ 
voir gas supply enabled by stellar feedback at earlier times, not 
through major mergers. With a duty cycle of nearly a gigayear, 
our simulations show that the submillimetre-luminous phase of 
high-z galaxies is a drawn-out one that is associated with signifi¬ 
cant mass buildup in early Universe proto-clusters, and that many 
submillimetre-luminous galaxies are actually composed of numer¬ 
ous unresolved components (for which there is some observational 
evidenc^^. 

We have conducted our cosmological hydrodynamic zoom galaxy 
formation simulations utilising the new hydrodynamic code GlZMCpl 
and include a model for the impact of radiative and thermal pressure 
from stars on the multiphase interstellar medium. This feedback both 
regulates the star formation rate, and shapes the structure in the in¬ 
terstellar medium. Informed by clustering measurements of observed 
smgP, we focus on a massive (Mdm ~ 10 ^^ M© at z = 2) halo 
with baryonic particle mass Mbary ~ 10^ M© as the host of our “main 
galaxy”, and run the simulation to = 2. The only condition of the 
tracked galaxy pre-selected to match the physical properties of ob¬ 
served SMGs is the chosen halo mass. We combine this with a new 
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dust radiation transport package, POWDERDAY, that simulates the tra¬ 
verse of stellar photons through the dusty ISM of the galaxy, allowing 
us to robustly translate our hydrodynamic simulation into observable 
measures. We simulate the radiative transfer from a 200 kpc region 
around the main galaxy. This simulation represents the first cosmo¬ 
logical model of a galaxy this massive to be explicitly coupled with 
dust radiative transfer calculations. The details of both codes and the 
simulation setup are fully described in the Methods section. 

We define two distinct regions in the simulations. The “submil¬ 
limetre emission region” is the 200 kpc region surrounding the cen¬ 
tral galaxy in the halo of interest. This is the region where all of the 
modeled 850/im emission comes from, and is what relates most di¬ 
rectly to observations. The “submillimetre galaxy” refers to the central 
galaxy in the halo. Physical quantities from the submillimetre galaxy 
are most applicable to high-resolution observations, as well as placing 
these models in the context of other theoretical galaxy formation mod¬ 
els. As we will show, the submillimetre emission from the region is 
generally dominated by the central submillimetre galaxy, though the 
contribution from lower mass galaxies is often non-negligible. 

We track the submillimetre properties of the galaxies within the re¬ 
gion from z 6. The star formation rates of galaxies in the region 
rise from this redshift toward later times z ^ 2, owing to accretion of 
gas from the intergalactic medium (Figure 1). As stars form, stellar 
feedback-driven galactic winds generate outfiows and fountains allow¬ 
ing recycled gas to be available for star formation at later times (ED Fig 
1). This phenomena shapes a star formation history that is still rising at 
z ~ 2, in contrast to galaxy formation models with more traditional im¬ 
plementations of su bresol ution feedback, which peak at z ~ 3 — 6 for 
galaxies of this mas^^®^ Mergers and global instabilities drive short¬ 
term variability in the global SFR, while outfiows and infall driven by 
the feedback model can impact features in the star formation history 
with a somewhat cyclical ’saw-tooth’ pattern. 

At its earliest stages (z ~ 4—6), the integrated SFR from the galax¬ 
ies in the region varies from ~ 100 — 300 M© yr“^, with a significant 
stellar mass (0.5 — lx 10^^M©) in place, comparable to some high- 
redshift detectionFeedback from massive stars enrich the interstel¬ 
lar medium with metals, and the dust content simultaneously rises. By 
z ^ 3, the combination of gas accumulation with substantial metal 
enrichment drives a factor of ~ 50 increase in the dust mass, with 
masses approaching ~ 1 x lO^M©. Radiation from the delayed peak 
in the star formation rate interacting with this substantive dust reser¬ 
voir drives the observed 850/im fiux density to detectable values of 
> 5 mJy. The galaxies associated with the main halo enter a long-lived 
submillimetre-luminous phase, with a duty cycle of ~ 0.75 Gyr. While 
our main model is only run to z = 2 owing to computational restric¬ 
tions for models of this resolution, tests with lower-resolution models 
reveal that at later times (z < 1.5), a declining star formation rate due 
to inefficient accretion as well as exhausted gas supply drives a drop 
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in the submillimetre flux density (for more details, see the Methods 
section). The star formation history of galaxies residing in (z = 2) 
Mdm ~ 10^^ Mq halos as controlled by the underlying stellar feed¬ 
back provides a physical explanation for the peak in the observed SMG 
redshift distribution at z = 2 — 

During the submillimetre-luminous phase, the emitting region is 
almost always occupied by multiple detectable galaxies. In Figure 2, 
we present gas surface density projections of six arbitrarily chosen 
snapshots during the evolution of the submillimetre-luminous phase 
(z = 2 — 3). The panels are 250 kpc on a side; for reference, the 
full width at half maximum (FWHM) of the Submillimetre Common- 
user Bolometer Array (SCUBA) on the James Clerk Maxwell Tele¬ 
scope (JCMT), the first instrument to detect SMGs, is ~ 125 kpc at 
z ^ 2. Multiple clumps of gas falling into the central are nearly always 
present. The observed flux density from the region is typically domi¬ 
nated by the central, with (on average) ~ 30% arising from emission 
from subhalos (ED Fig 2). The submillimetre flux density of the central 
galaxy rises dramatically between z ^ 2 — 3, reaching a peak value of 
~ 20 mJy. Owing to contributions from subhalos surrounding the cen¬ 
tral, the flux from the overall 200 kpc region can exceed this, peaking 
at ~ 30 mJy. Similarly extreme systems have recently been detecte d 
with the Herschel Space Observatory and South Pole Telescop^EMll 

While the central is being bombarded by subhalos over a range 
of mass ratios during the submillimetre-luminous phase, major galaxy 
mergers akin to local prototypical analogues such as Arp 220 or NGC 
6240 do not drive the onset of the long-lived submillimetre-luminous 
phase in the central galaxy. In Figure 1, we highlight when the galaxy 
undergoes a major merger with mass ratio >1:3. While major merg¬ 
ers are common at early times (and indeed drive some short-lived bursts 
in star formation), the bulk of the submillimetre-luminous phase at later 
times (z ^ 2 — 3) occurs nearly a gigayear after the last major merger. 
The ratio of the SFR to its integral over cosmic time (the specific SFR) 
of the overall emitting region is generally on the main sequence of 
galaxy formation at z ~ 2 (defined as the main locus of points on 
the SFR-Mhc relation), though the central can have values comparable 
both to main sequence galaxies between z ~ 2 — 3, as well as out¬ 
liers. One consequence of a model in which SMGs typically lie on the 
main sequence of star formation is that the gas surface densities show 
a broad range, from ~ 100 — lO^Mopc”^ (ED Eig 3), as well as di¬ 
verse gas spatial extents (Eigurej^. This is manifested observationally 
in the broad swath occupied by SMGs in the Kennicutt-Schmidt star 
formation relatioiPl The gaseous spatial extent and surface densities is 
to be contrasted, however, with local merger-driven ultraluminous in¬ 
frared galaxies (ULIRGs), which exhibit typical full width at half max¬ 
imum radii of ~ 100 — 500 pcP^. Idealised galaxy merger simulations 
with initial conditions designed to form SMGs further underscore this 
contrast, as they also result in compact morphologies during final co¬ 
alescence, and can be inefficient producers of submillimetre radiation 
owing to increased dust temperature^. 

The central submillimetre galaxy is amongst the most massive and 
highly star-forming of galaxies at this epoch. The stellar masses are 
diverse, and range from ~1 — 5x10^^ M©, comparable to recent 
measurements of this populatiorP^ as well as constraints from abun¬ 
dance matching techniquePl The molecular gas fractions of the cen¬ 
tral galaxy (/gas = Mh 2 /(Mh 2 + M*)) decline with stellar mass, and 
range from ~ 40% at lower stellar masses to < 10% at the highest 
masses. This is in agreement with observation^, though is depen¬ 
dent on the conversion from carbon monoxide (^^CO) luminosity to 
H 2 gas mass. We note that these predictions are quantitatively differ¬ 
ent from those produced by previous cosmological efforts in this held, 
with some predicted gas fractions exceeding /gas = 0.7^0 

and me¬ 
dian stellar masses as low as ~ 10^° M©^. We present the plots high¬ 
lighting the gas fractions, and calculated SEDs of our model SMG in 


the context of observations in ED Eigs 4-5. The gas distributions within 
the central galaxy, which range from ~ 1 — 8 kpc, compare well with 
recent observed dust maps with the Atacama Large Millimetre Arra}P^. 

The stellar masses, gas fractions and duty cycles are in agreement 
with previous lower-resolution cosmological effortJ^, though the pre¬ 
dicted SER and luminosity from this model are substantially larger. 
The star formation rate of the group of galaxies in the region peaks at 
~ I 5 OOM 0 yr“^. Importantly, up to half of the total infrared lumi¬ 
nosity can come from older stars with ages tage > O.lGyr. Utilising 
standard conversionP^, the estimated star formation rate from the in¬ 
tegrated infrared SED (3 — llOO/xm) can exceed ~ 3000Mo yr“^ 
(ED Eig 6), and hence infrared-based star formation rate derivations of 
dusty galaxies at high-z may over-estimate the true SER by a factor 
~ 2. Indeed, the contribution of satellite galaxies to the global SER, 
alongside the contribution of old stars to the infrared luminosity may 
relieve some tensions between the inferred star formation rates from 
submillimetre galaxies and massive galaxies modeled in cosmological 
hydrodynamic simulations^. 

The end-product of the central submillimetre galaxy at 2 : ~ 2 is a 
galaxy with a stellar mass of ~ 4 — 5 x 10^^ M©, distributed over a 
similarly compact region of ~ 1 — 5 kpc as the gas (Eigurej^. This is 
similar in extent and mass to the z ~ 2 compact quiescent galaxy popu¬ 
lation, an observed population with mean half-light radius of ~ 1-5 
kpc, stellar mass M* > 10^and ages tage ~ 0.5 — IGyP, 
suggesting a plausible connection between the galaxy populations. In¬ 
deed, a calculation of the stellar velocity dispersion along three orthog¬ 
onal sightlines of the central during the submillimetre-luminous phase 
results in cr* 600 — 700 km s“^, comparable to measurements of 
high-z compact quiescents. A large sample of simulated SMGs will al¬ 
low for a robust analysis of expected abundances of SMGs and compact 
quiescents at z ~ 2. 

Our picture for SMG formation suggests that they are not transient 
events, but rather natural long-lived phases in the evolution of mas¬ 
sive halos. The ~ 0.75 Gyr duty cycle combined with the comoving 
abundancJ^ of dark matter halos of this mass results in an expected 
abundance of our model SMGs of~ 1.5 x 10~^h^ Mpc“^, compara¬ 
ble to the ~ Mpc“^ observed for SMgJ^. While modeling the 

full number counts involves convolving the typical duty cycle as a func¬ 
tion of halo mass with halo mass functions over a range of redshifts, the 
approximate abundances implied by this model are encouraging. 

This model suggests that galaxies that form in halos of mass 
Mdm ~ 10^^ at z = 0 will represent typical SMGs near the peak 
of their redshift distribution. Lower mass models show that they do not 
achieve the requisite star formation rate and metal enrichment to gen¬ 
erate submillimetre-luminous galaxies (see Metho ds sec tion). More 
extreme SMGs being detected between z = 5 — may form in 
even more massive (and rare) halos than those considered here. 


1. Casey, C. M., Narayanan, D. & Cooray, A. Dusty star-forming galaxies at 
high redshift. Physics Reports 541 , 45-161 (2014). 

2. Hickox, R. C. et al. The LABOCA survey of the Extended Chandra Deep 
Field-South: clustering of submillimetre galaxies. Mon. Not. R. Astron. 
Soc. 421,284-295 (2012). 

3. Geach, J. E. et al. On the Evolution of the Molecular Gas Fraction of 
Star-Forming Galaxies. Astrophys. J. Let. 730 , LI 9 (2011). 

4. Fu, FI. etal. The rapid assembly of an elliptical galaxy of 400 billion solar 
masses at a redshift of 2.3. Nature 498 , 338-341 (2013). 

5. Daddi, E. etal. Two Bright Submillimeter Galaxies in a z = 4.05 Protoclus¬ 
ter in Goods-North, and Accurate Radio-Infrared Photometric Redshifts. 
Astrophys. J. 694 , 1517-1538 (2009). 

6. Swinbank, A. M. etal. The Rest-Frame Optical Spectra of SCUBA Galax¬ 
ies. Astrophys. J. 617 , 64-80 (2004). 

7. Baugh, C. M. etal. Can the faint submillimetre galaxies be explained in 
the A cold dark matter model? Mon. Not. R. Astron. Soc. 356 , 1191- 
1200 (2005). 


2 




8. Hayward, C. C. et al. Submillimetre galaxies in a hierarchical universe: 
number counts, redshift distribution and implications for the IMF. Mon. 
Not. R. Astron. Soc. 428, 2529-2547 (2013). 

9. Shimizu, I., Yoshida, N. & Okamoto, T. Submillimetre galaxies in cosmo¬ 
logical hydrodynamic simulations: source number counts and the spatial 
clustering. Mon. Not. R. Astron. Soc. 427, 2866-2875 (2012). 

10. Dave, R. et al. The nature of submillimetre galaxies in cosmological 
hydrodynamic simulations. Mon. Not. R. Astron. Soc. 404, 1355-1368 
( 2010 ). 

11. Hodge, J. A. et al. An ALMA Survey of Submillimeter Galaxies in the 
Extended Chandra Deep Field South: Source Catalog and Multiplicity. 
Astrophys. J. 768, 91 (2013). 

12. Hopkins, R F. GIZMO: A New Class of Accurate, Mesh-Free Hydrody¬ 
namic Simulation Methods. arXiv/1409.7395 (2014). 

13. Dave, R., Finlator, K. & Oppenheimer, B. D. An analytic model for the 
evolution of the stellar, gas and metal content of galaxies. Mon. Not. R. 
Astron. Soc. 421, 98-107 (2012). 

14. Hopkins, P. F. et al. Galaxies on FIRE (Feedback In Realistic Environ¬ 
ments): stellar feedback explains cosmologically inefficient star forma¬ 
tion. Mon. Not. R. Astron. Soc. 445, 581-603 (2014). 

15. Feldmann, R. & Mayer, L. The Argo Simulation: I. Quenching of Mas¬ 
sive Galaxies at High Redshift as a Result of Cosmological Starvation. 
arXiv/1404.3212 (2014). 

16. Finkelstein, S. L. et al. A galaxy rapidly forming stars 700 million years 
after the Big Bang at redshift 7.51. Nature 502, 524-527 (2013). 

17. WeiB, A. et al. ALMA Redshifts of Millimeter-selected Galaxies from the 
SPT Survey: The Redshift Distribution of Dusty Star-forming Galaxies. 
Astrophys. J. 767, 88 (2013). 

18. Ivison, R. J. etal. Herschel-ATLAS: A Binary HyLIRG Pinpointing a Clus¬ 
ter of Starbursting Protoellipticals. Astrophys. J. 772, 137 (2013). 

19. Hezaveh, Y. D. et al. ALMA Observations of SPT-discovered, Strongly 
Lensed, Dusty, Star-forming Galaxies. Astrophys. J. 767, 132 (2013). 

20. Downes, D. & Solomon, P M. Rotating Nuclear Rings and Extreme Star- 
bursts in Ultraluminous Galaxies. Astrophys. J. 507, 615-654 (1998). 

21. Michafowski, M. J. et al. The stellar masses and specific star-formation 
rates of submillimetre galaxies. Astron. Astrophys. 541, A85 (2012). 

22. Behroozi, P S., Wechsler, R. H. & Conroy, C. The Average Star 
Formation Histories of Galaxies in Dark Matter Halos from z=0-8. 
arXiv/1207.6105 (2012). 

23. Tacconi, L. J. etal. Phibss: Molecular Gas Content and Scaling Relations 
in z ~ 1-3 Massive, Main-sequence Star-forming Galaxies. Astrophys. J. 
768, 74 (2013). 

24. Simpson, J. M. et al. The SCLIBA-2 Cosmology Legacy Survey: ALMA 
Resolves the Rest-frame Far-infrared Emission of Sub-millimeter Galax¬ 
ies. Astrophys. J. 799, 81 (2015). 

25. Kennicutt, R. C. & Evans, N. J. Star Formation in the Milky Way and 
Nearby Galaxies. Annu. Rev. Astron. Astrophys. 50, 531-608 (2012). 

26. van Dokkum, P G. et al. Confirmation of the Remarkable Compactness 
of Massive Quiescent Galaxies at z ~ 2.3: Early-Type Galaxies Did not 
Form in a Simple Monolithic Collapse. Astrophys. J. Let. 677, L5-L8 
(2008). 

27. Murray, S. G., Power, C. & Robotham, A. S. G. HMFcalc: An online tool 
for calculating dark matter halo mass functions. Astronomy and Comput¬ 
ing 3, 23-34 (2013). 

28. Chapman, S. C., Blain, A. W., Smail, I. & Ivison, R. J. A Redshift Survey of 
the Submillimeter Galaxy Population. Astrophys. J. 622, 772-796 (2005). 

29. Riechers, D. A. et al. A dust-obscured massive maximum-starburst 
galaxy at a redshift of 6.34. Nature 496, 329-333 (2013). 

30. Vieira, J. D. et al. Dusty starburst galaxies in the early Universe as re¬ 
vealed by gravitational lensing. Nature 495 , 344-347 (2013). 

Acknowledgements The authours thank Michaf J. Michafowski for provid¬ 
ing observational data. Partial support for DN was provided by NSF AST- 
1009452, AST-1442650, NASA HST AR-13906.001, and a Cottrell College 
Science Award. PFH, CCH, MT and RT were funded by the Gordon and 
Betty Moore Foundation (GBMF4561 and Grant #776). PFH acknowledges 
the Alfred P Sloan Foundation for support. CAFG was supported by NASA 
awards PF3-140106, NNX15AB22G, and NSF AST-1412836. DK was sup¬ 
ported by NSF AST-1412153. RF was supported by NASA HF-51304.01-A . 
The simulations here were run on Stampede at TACC through NSF XSEDE 
allocations #TG- ASTI 20025, TG-AST130039, and TG-AST140023, NASA 
Pleiades, and the Haverford College cluster. 

Author Contributions D.N. wrote the text, and led the radiative transfer sim¬ 
ulations and analysis. D.N., M.T, TR. and R.T wrote the Powderday soft¬ 
ware. R.T, C.C.H. and D.B. contributed to simulation analysis, and R.F, P.H., 
C-A.F-G and D.K. performed the cosmological simulations 


Author Information Reprints and permissions information is available at 
www.nature.com/reprints. The authors declare no competing financial inter¬ 
ests. Readers are welcome to comment on the online version of the pa¬ 
per. Correspondence and requests for materials should be addressed to D.N. 
(dnarayan@haverford.edu). 


3 



Time since Big Bang (Gyr) Time since Big Bang (Gyr) 



Figure 1: Evolution of physical and observable properties of submillimetre emission region and central galaxy. The 200 kpc submillimetre 
emission region are shown with thick solid lines, while the central galaxy’s properties are given by thin dashed lines in each panel. Stellar and 
dust mass are in the top left, SFR at the top right; predicted observed 850 /xm flux density at the bottom left; and specific SFR (M*/SFR) at the 
bottom right. The SFR is averaged on 50 Myr timescales, and includes a correctional factor 0.7 for mass loss. Locations of major galaxy mergers 
(> 1 : 3) are noted by green vertical ticks on the top axis of the top right panel. The blue shaded region in the 850 /xm curve shows when the 
galaxy would be detectable as an SMG with SCUBA (iSsso > 5 mJy). The yellow and purple shaded regions in the bottom right show the rough 
ranges for the = 2 Main Sequence and Starburst regime. The grey region denotes below the Main Sequence. 
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Figure 2: Surface density projection maps of 250 kpc region around central submillimetre galaxy between redshifts 2 — 3. The 

submillimetre-emission region probed in surveys typically encompasses a central galaxy in a massive halo that is undergoing a protracted 
bombardment phase by numerous sub-halos. Some of the brightest SMGs arise from numerous galaxies within the beam in a rich environment 
(bottom right panel). 
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Central Galaxy Radius (kpc) 

Figure 3: Gas and stellar radius distribution for central submillimetre galaxy. The orange histogram denotes the half mass radius of the stars, 
while the blue shows the gas. The galaxy gas is more distributed in the central than the (sub-kiloparsec) extent expected from major mergers, 
though still sufficiently compact that it will remain unresolved even with ~ arcsecond resolution. 
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METHODS 


1 Cosmological Hydrodynamic Zoom Simulations 

We utilise a newly-developed version of TreeSPH that employs 
a pressure-entropy formulation of smoothed particle hydrodynamics 
(sphP that obviates many of the potential discrepancies noted be¬ 
tween gri d-bas ed codes, traditional SPH codes, and moving-mesh 
algorithmic^®. In particular, we employ the hydrodynamic code 
GiZMcP in P-SPH mode which conserves momentum, energy, angu¬ 
lar momentum and entropy, and includes newly developed algor ithms 
to treat the artificial viscosity, entropy diffusion and time-steppingp“ 
The gravity solver is a modified version of the GADGET-3 solvei“ and 
an updated softening kernel to better represent the potential of the SPH 
smoothing kernel is includecP^, 

The simulations are fully cosmological zoom-in calculations of the 
evolution of individual galaxies. A 144 Mpc^ cosmological volume 
was simulated at low resolution down to redshift z = 0 with dark mat¬ 
ter only. The halo of interest was identified, and re-simulated at much 
higher resolution with baryons included. The initial conditions were 
generated with the MUSIC cod^. We simulate four zoom galaxies - 
one is our main galaxy, and the other three are at varying resolutions 
and masses for the purposes of testing. The main galaxy of interest to 
this study resides in a dark matter halo mass of Mdm = 3 x 10^^ M© 
at z = 2. The initial baryonic particle masses in the high-resolution 
region were 2.7 x 10^ Mq, and the minimum baryonic/stars/dark mat¬ 
ter force softening lengths were 9/21/142 proper pc at z = 2. The 
physical properties of all of the modeled galaxies are presented in ED 
Table 1 in the Extended Data. 

The baryonic physics implemented into GIZMO are developed 
based on extensive tests studyi ng idealised simulations of both isolated 
disks and galaxy mer^rs^^®. The gas cools utilising an updated cool¬ 
ing curve to standarcP^ implementations in SPH codes which includes 
both atomic and molecular line emissiorl^. The modeled interstel¬ 
lar medium is multiphase. The neutral ISM is broken into an atomic 
and molecular component following algorithms that scale the m olec- 
ular fraction with column density and gas phase metallicit}^®. Star 
formation occurs in molecular gas above a threshold density (here, this 
is set to nthresh = 10 cm“^). Star formation is further restricted to gas 
that is locally self-gravitating, where: 


a = ^' 


IV-V 


" + |Vx v|" 


Gp 


< 1 


( 1 ) 


This follows from studieP^ that show that the predicted spatial distri¬ 
bution of star formation in galaxies is more realistic when utilising a 
gas self-gravitating criterion compared to a variety of other algorithms 
(including a fixed density threshold, a pure molecular-gas law, a tem¬ 
perature threshold, a Jeans criterion, a cooling-time criterion and a con¬ 
verging flow criterion). The star formation rate follows a volumetric 
relation: 

P* = Pmol / (2) 

That is, stars are allowed to form with 100% efficiency per free fall 
time. The star formation is subsequently self-regulated by stellar feed¬ 
back, resulting in a time-averaged efficiency on galaxy scales of eff of 
~ 0.005 - O.ll^. 

Once stars have formed, they can impact the ISM via various feed¬ 
back mechanisms. Assumi ng a Kroup#^ stellar initial mass function, 
and utilising STARBURST99^ for luminosity, mass-return and super¬ 
nova rate calculations as a function of stellar age and metallicity, we 
include the following forms of stellar feedback: 

1. Radiation Momentum Deposition: At each timestep, the gas near 
young stars is impacted by a momentum fiux given by 

Prad ~ exp( '7"uv/optical)) (1 T "^IR)-^^incident/C (3) 


Where tir is calculated directly from the simulation, as tir = 
Sgas/^iR, and a^ir = 5(Z/Zo) g“^ cm^. 

2. Supernovae and Stellar Winds : We utilise tabulated Type-1 and 
Type-II supemovae rateJ2E3 if a supernova occurs during a 
timestep, thermal energy and radial momentum are injected within 
a smoothing length of the star. Gas and metal return is included as 
well. Stellar winds are similarly included with energy, wind mo¬ 
mentum, mass and metals deposited within a smoothing length. 

3. Photoheating of HII regions: The production rate of ionising ra¬ 
diation from stars determines the extent of HII regions (allowing 
for overlapping regions). These regions are heated to 10^ K if the 
gas is below that threshold. 

We utilise the second and third models in the Table to test the con¬ 
vergence properties of our simulations. One model is run with the same 
mass baryonic resolution as our main model (standard resolution; SR), 
and one a factor of ~ 8 higher resolution (high resolution; HR). In ED 
Eig 7, we show the modeled duty cycle above a given flux density as a 
function of flux density for these two models. We see that the shortest 
lived (< 200 Myr) emission spikes present in the standard resolution 
model may not be converged in the highest resolution model. Notably, 
emission with longer duty cycles are either converged, or underpre¬ 
dicted in our standard resolution model, suggesting that the relatively 
long-lived submillimetre-luminous phase is robust. 

We show the M* — z relation for the central in ED Eig 8 as com¬ 
pared to observational constraints^. The central galaxy has a stellar 
mass a factor ~ 2 greater than the observed median stellar mass for 
comparable mass halos at this epoch. The model galaxy may represent 
an outlier in the M* — z relation at these redshifts. Indeed, the thick¬ 
ness of the observational constraints shows the uncertainty, not range in 
possible values. Alternatively, it is possible that the inclusion of feed¬ 
back from an active galactic nucleus (AGN) could impact the stellar 
mass buildup in the galaxy, though the level to which black hole growth 
can impact star formation near the submillimetre-luminous phase is un¬ 
clear. Some models have show n that AGN can grow efficiently in the 
absence of major mergerJSEI while other models and ob serva tions 
suggest that mergers may be necessary to grow massive holeP^®. The 
last major merger before the submillimetre-luminous phase is ~ 1 Gyr 
prior. Tests with our low resolution model (ml3ml4) show that with¬ 
out AGN feedback, residual star formation drives a factor ~ 2 increase 
in stellar mass at late times (z 0 — 1). Einally, we note that a higher 
mass resolution model could potentially also result in decreased final 
stellar masses. In our convergence tests, the final (at z — 2) of the 
HR run is ~ 60% that of the SR run. 

2 Dust Radiative Transfer Caiculations 

To calculate the inferred observational properties of our simulated 
galaxies, we have developed a dust radiative transfer package, POW- 
DERDAY. In short, POWDERDAY takes hydrodynamic simulations of 
galaxies in evolution, projects the gas properties onto an adaptive mesh, 
and calculates the radiative transfer from the stellar sources through 
the dusty interstellar medium until an equilibrium dust temperature is 
achieved. 

In detail, we identify galaxi es utilising SKID to locate bound 
gro ups o f baryonic particleand track their progenitors back in 
tima^^I^. Galaxies and halos are required to contain at least 64 par¬ 
ticles each in order to be identified. We cut out a 200 kpc (side 
length) region around the galaxy of interest, and subdivide the do¬ 
main into an adaptive grid with an octree memory structure. Eor- 
mally, we begin with one cell encompassing the entire 8 x 10® kpc^ 
radiative transfer region. The cells then recursively subdivide into 
octs until there are a threshold maximum number of gas particles in 
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the cell (we employ risubdivide,thresh = 64, though experiments with 
^subdivide,thresh = 32 show convcrged results). The physical prop¬ 
erties of the gas particles are projected onto the octree using a spline 
smoothing kerneP- 

The spectral energy distribution of stars are calculated o n the 
fly with the Flexible Stellar Population Synthesis code, 
through PYTHON-FSPS, a set of PYTHON hooks for FSPS (https : 
//github . com/dfm/python-fsps ). The SEDs are calculated 
as simple stellar populations with ages and metallicities determined by 
the hydrodynamic simulation, and assuming a Kroupa IMF. 

The radiative transfer happens in a Monte Carlo fashion utilising 
the three-dimensional dust radiative transfer solver, HYPERlOfJ^. The 
code uses an iterative methodology to determine the radiative equilib¬ 
rium temperaturJ^, and we determine convergence when the energy 
absorbed by 99% of the cells has changed by less than 1% between it¬ 
erations. We assume a dust grain-size distribution comparable to that of 
the Milky Wa 5 l^ with R = AvlE{B — y) = 3.15. The dust emissiv- 
ities are updated to include an approximation for polycyclic aromatic 
hydrocarbons (PAHs) alongside thermal emission^. We assume a con¬ 
stant dust to metals ratio of 0.4, m otivate d by both Milky Way and 
extragalactic observational constraintP^^. 

The underlying HYPERION code has passed the standard bench¬ 
marks for codes of this typJ^, and we have found that POWDERDAY 
com pares well against other publicly available dust radiative transfer 
codeP^^^in test starburst SPH galaxy merger simulations. 

3 Parameter Choices 

In ED Fig 9, we present a number of tests of our parameter choices 
for the radiative transfer calculations. We show the predicted 850/im 
light curve from our lowest resolution model (ml3m 14) utilising both 
flducial parameters, as well as three parameter choice variations. 

We first ask whether our chosen radiative transfer grid size affects 
our principle results. Our flducial model is a 200 kpc (on a side) box 
cut out of the global cosmological simulation centred on the halo of 
interest. This size was chosen to reflect a rough average of the typical 
(sub)millimetre beam sizes typically used to detect SMGs. For ex¬ 
ample, assuming Planck 2013 cosmological parameterP^, the Submil¬ 
limetre Common-Use Bolometer Array (SCUBA) on the James Clerk 
Maxwell Telescope (JCMT) has a 15 full width at half maximum 
(FWHM) beam at 850/im. At z = 2 this corresponds to ~ 128 
kpc. At the same redshift, the beam of AzTEC and LAB OCA at 1 
mm on the JCMT corresponds to ~ 163 kpc (19‘‘); The South Pole 
Telescope (SPT) has a beamsize of 540 kpc at 1.4 mm (63 ); and Her- 
schel’s SPIRE instrument ranges from 154 — 308 kpc (250 — 500 /xm; 
18 - 36“). 

Because a few notable beam sizes (of particular relevance, the 
SCUBA beam) are smaller than our assumed box size of 200 kpc, we 
have run an additional model with box length 100 kpc (all other pa¬ 
rameters exactly the same). We highlight the resultant 850 /xm light 
curve from this model in the top right panel of ED Fig 9. When com¬ 
paring to our flducial model, it is apparent that our results are robust to 
the highest resolution beams that have been used for SMG surveys at 
single dish facilities to date. 

We additionally investigate whether our inclusion of PAHs in our 
model makes any difference to the calculated submillimetre-wave flux 
density of our model galaxy. This is presented in the bottom left panel 
of ED Fig 9. Again, we note minimal impact on the submillimetre SED 
of our model. 

Finally, we ensure that our results are converged with the number 
of photons emitted. We flducially run 10^ photons per grid (roughly 
100 per cell). In the bottom-right, we show the results from a run with 
10® photons per grid, and show that the results are robust against this 
parameter choice. 


4 Relation to other Models 

Historically, the methods used, and physical models for SMG for¬ 
mation in numerical simulations are quite varied. Here, we summarise 
these methods and results, and place our own model into this context. 
Broadly, there are three classes of SMG formation models: cosmologi¬ 
cal semi-analytic models (SAMs), idealised non-cosmological simula¬ 
tions, and cosmological hydrodynamic models. The present model falls 
into the latter category. Our presented model is the first self-consistent 
cosmological simulation with baryons and bona fide radiative transfer 
to form a submillimetre galaxy with physical properties comparable to 
those observed. 

The initial forays into this held were typically with SAMs. This 
owes to the fact that SAMs are computationally inexpensive, and allow 
for a large search in physical parameter space relatively easily. SAMs 
either utilise analytic halo merger trees, or directly simulate them, and 
then employ a naly tic prescriptions to describe the central galaxies. The 
Durham SAIVP^ couples galaxies formed in a semi-analytic model 
with dust radiative transfer. These simulations model galaxies that have 
axisymmetric geometries that consist of a disc and a bulge. Young stel¬ 
lar populations are assumed to still be enveloped in their birth clouds, 
and thus experience additional attenuation. This model suggests that 
SMGs typically owe to ~ 22% major mergers, the remainder as mi¬ 
nor mergers, and that the stellar initial mass function is flat during the 
starburst. The typical duty cycle for the submillimetre-luminous phase 
is ~ 100 Myr (a factor ~ 7.5 lower than found in our work), galaxies 
extremely gas rich (/gas ~ 75%), and stellar masses a factor ~ 10 
lower than our mode l (M^ ^ 2. x lO^^M©). While the stellar masses 
of SMGs are debate J^ ^^l^^ j the gas fractions appear to be uniformly 
lower in observationJ3i3S®Sl, and a flat stellar initial mass function 
likely ruled out by CO dynamical mass measurementJ^. 

As an alternative to cosmological simulations, a nu mber o f works 
have explored SMG formation in idealised simulationJSSEl These 
works evolve hydrodynamic models of idealised discs and mergers 
over a range of merger mass ratios, and combine these with dust radia¬ 
tive transfer simulation^. These models infer halo masses and stellar 
masses for SMGs comparable to those modeled here. This said, in 
the idealised galaxy models, ~ 30 — 70% of the SMGs (flux depen¬ 
dent) owe to merger-driven starbursts, substantially higher than what 
is found for our model. Some workP^ have noted that binary mergers 
that cause SMGs may break up into multiples at high-resolution owing 
to the contribution to the total flux of individual inspiralling discs. Be¬ 
cause idealised simulations are non-cosmological in nature, comparing 
the multiplicity inferred from these to our models is difficult: the ma¬ 
jor merger multiplicity can only be 2 when considering galaxies at the 
same redshift. On the other hand ED Fig 2 suggests potentially larger 
multiplicity can be observed for physically associated clumps. 

To fully capture the cosmic environment of SMGs in formation, 
as well as their baryonic structure and morphology, cosmological hy¬ 
drodynamic simulations are likely the best tool. Thusfar, cosmological 
hydrodynamic simulations used to simulate SMGs have not employed 
direct radiative transfer modelJ^. As such, inferring when a galaxy is 
an SMG in cosmological simulations has necessitated the use of param- 
eterised emission models, such as assumed grey-body emission lawP, 
or star formation rate thresholdJ^. The physical properties for SMGs 
derived from the most extensive of these workJ^(i.e. M*, Mdm, and 
/gas) are similar to the model presented here, though with roughly a 
factor ~ 3 difference in SFR. 

5 Code Availability 

We have made POWDERDAY available at https: 
//bitbucket. org/desika/powderday . and GIZMO 
available at https : //bitbucket. org/phopkins/gizmo 
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Table Extended Data 1: Summary of Model Galaxies. M* and Mhaio refer to the stellar and halo mass at z — 2] Cb and cbm refer to the 
minimum force softening lengths for baryons and dark matter particles. 
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Extended Data Figure 1: Mass of gas in central galaxy that will be consumed during SMG starburst as a function of z. The gas mass 
consumed during the starburst is calculated by tracking the evolution of gas particles that turn into stars during the SMG phase (z ^ 2 — 2.7), 
and is only measured for the central galaxy itself (i.e. gas ejected into the halo is not included). The colours denote the median scale height from 
the galaxy centre of mass. The SMG gas reservoir follows a cycle of being pushed outward followed by re-accretion. 
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Extended Data Figure 2: Predicted contribution of submillimetre-luminous region components to total flux density. Submillimetre- 
luminous regions often break up into multiples. Shown is a histogram of the ratio of the brightest component to the total flux density from 
the region, with the average denoted with the vertical line. The region is generally dominated by one component, though smaller subhalos can 
contribute on average ~ 30% of the observed flux density. 
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Extended Data Figure 3: Gas surface density for central submillimetre galaxy. The blue histograms show the time-weighted distribution of 
surface densities during all phases, whilst the orange show the same for the submillimetre-luminous phase. We predict that the submillimetre- 
luminous phases do not have drai uatic ally different surface density distributions compared to the non-submillimetre-luminous phases. This 
feature may be tentatively observecl®!. 
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Extended Data Figure 4: Molecular gas fraction as a function of galaxy stellar mass. Blue stars show individual snapshots of central 
submillimetre galaxy whilst red circles with error bars (la) show observations with direct CO (J=l-0) measurements (to avoid complications in 
converting from higher-lying CO rotational lines to the ground state for a mass conversion). Both observations and our model show a declining 
molecular gas fraction with increasing galaxy mass, with a typical range of /gas = 0.1 — 0.4 for galaxies of SMG mass. 
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Extended Data Figure 5: Predicted spectral energy distribution (SED) for central submillimetre galaxy. The blue shaded region shows the 
range of SEDs for all snapshots that satisfy the fiducial FssoAxm > 5 mJy submillimetre galaxy selection criteria, whilst the grey points with error 
bars (la) are a compilation of observed data. The individual coloured lines show the SEDs for individual submillimetre luminous snapshots. 
The data and models are redshifted to a common redshift z = 2. The model and data compare well, and the model suggests a diverse range of 
SMG SEDs. 


5 




























SFR5o(Mgyr-M 

Extended Data Figure 6: Overestimate of SFR of High-z SMGs. The ordinate denotes the SFR as determined from the infrared seeP, while 
the abscissa shows the SFR averaged over the last 50Myr in the simulations. Up to an SFR of ~ 8 OOM 0 yr“^ the two correspond well. At higher 
SFRs, however, there is a dramatic departure owing to substantial contribution to the infrared luminosity by older stars. 
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Extended Data Figure 7: Resolution tests for hydrodynamic zoom simulations. Lines show 850/im duty cycle above a given flux density as 
a function of flux density for our resolution test models presented in the Methods section. SR denotes our standard resolution (the resolution of 
our main model) whilst HR is one level higher reflnement. 
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Extended Data Figure 8: Stellar Mass - Redshift Relation for Model Galaxy. Purple line shows model whilst blue filled region shows 
observational constraints from an abundance matching assumptioiP^. The model and observations are in reasonable agreement, especially during 
the submillimetre-luminous phase (purple shaded region). At late times, the stellar mass of the galaxy is a factor ~ 2 higher than the median 
observed galaxy. 
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Extended Data Figure 9: Tests of parameter choices for radiative transfer calculations. The simulated galaxy for these tests is our lowest 
resolution cosmological simulation (ml3ml4). Each panel shows the 850 fim flux density lightcurve of the tested model, with time noted on the 
abscissa (redshift on the bottom, time since the Big Bang on top). In all panels, the purple shaded region denotes S'sso > 5 mJy, the canonical 
selection criteria for SMGs. Top Left: Our flducial set of parameters; Top Right: Simulation with a 100 kpc (on a side) emission region instead 
of 200 kpc; Bottom Left: Simulation with our model for PAHs turned off; Bottom Right: Fiducial simulation run with ten times the number of 
photons. 
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